Calculation of material properties and ray tracing in transformation media 
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Complex and interesting electromagnetic behavior can be found in spaces with non-flat topology. 
When considering the properties of an electromagnetic medium under an arbitrary coordinate trans- 
formation an alternative interpretation presents itself. The transformed material property tensors 
may be interpreted as a different set of material properties in a flat, Cartesian space. We describe 
the calculation of these material properties for coordinate transformations that describe spaces with 
spherical or cylindrical holes in them. The resulting material properties can then implement invisi- 
bility cloaks in flat space. We also describe a method for performing geometric ray tracing in these 
materials which are both inhomogeneous and anisotropic in their electric permittivity and magnetic 
permeability. 



I. INTRODUCTION 
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Recently the use of coordinate transformations to pro- 
duce material specifications that control electromagnetic 
fields in interesting and useful ways has been discussed. 
We have described such a method in which the transfor- 
mation properties of Maxwell's equations and the consti- 
tutive relations can yield material descriptions that im- 
plement surprising functionality, such as invisibility [1]. 
Another author described a similar method where the 
two dimensional Helmholtz equation is transformed to 
produce similar effects in the geometric limit [2]. We note 
that these theoretical design methods are of more than 
academic interest, as the material specifications can be 
implemented with metamaterial technology [3-8]. There 
has also been recent work on invisibility cloaking that 
does not employ the transformation method[9, 10]. 

In this article we describe how to calculate these ma- 
terial properties directly as Cartesian tensors and to per- 
form ray tracing on devices composed of these materials. 
We work out two examples in detail, the spherical and 
cylindrical invisibility cloak. 



II. MATERIAL PROPERTIES 

We will describe the transformation properties of the 
electromagnetic, material property tensors. We first note 
that the Minkowski form of Maxwell's equations [11-13] 
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is form invariant for general space-time transformations. 
F a p is the tensor of electric field and magnetic induction, 
and G a ° is the tensor density of electric displacement and 
magnetic field, and is the source vecto. In SI units 
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All of the information regarding the topology of the space 
is contained in the constitutive relations 



1 



(3) 



where C a ^^ v is the constitutive tensor representing the 
properties of the medium, including its permittivity, per- 
meability and bianisotropic properties. C a/3flL ' is a tensor 
density of weight +1, so it transforms as[ll] 
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written in terms of the Jacobian transformation matrix 

dx a ' 
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which is just the derivative of the transformed coordi- 
nates with respect to the original coordinates. If we re- 
strict ourselves to transformations that are time invari- 
ant, the permittivity and permeability are also tensors 
individually. Specifically, they are tensor densities of 
weight +1, which transform as[ll, 14] 
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where the roman indices run from 1 to 3, for the three 
spatial coordinates, as is standard practice. Equations 
(6) are the primary tools for the transformation design 
method when the base medium is not bianisotropic and 
the desired device moves or change shape with speeds 
much less than that of light, i.e. probably all devices of 
practical interest. These equations can be shown to be 
exactly equivalent to the results derived by Ward and 
Pendry[15]. 

If the original medium is isotropic, (6) can also be writ- 
ten in terms of the metric[14, 16] 
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Maxwell's equations, (1), together with the medium 
specified by (4) or (6) describe a single electromagnetic 
behavior, but this behavior can be interpreted in two 
ways. 

One way, the traditional way, is that the material prop- 
erty tensors that appear on the left and right hand sides 
of (4) or (6) represent the same material properties, but 
in different spaces. The components in the transformed 
space arc different form those in the original space, due 
to the topology of the transformation. We will refer to 
this view as the topological interpretation. 

An alternative interpretation, is that the material 
property tensors on the left and right hand sides of (4) 
or (6) represent different material properties. Both sets 
of tensor components are interpreted as components in a 
flat, Cartesian space. The form invariance of Maxwell's 
equations insures that both interpretations lead to the 
same electromagnetic behavior. We will refer to this view 
as the materials interpretation. 

To design something of interest, one imagines a space 
with some desired property, a hole for example. Then 
one constructs the coordinate transformation of the space 
with this desired property. Using (4) or (6) one can then 
calculate a set of material properties that will implement 
this interesting property of the imagined space in our own 
boring, flat, Cartesian space. 



A. Spherical Cloak 

The spherical cloak is designed by considering a spher- 
ically symmetric coordinate transformation. This trans- 
formation compresses all the space in a volume of radius, 
b, into a spherical shell of inner radius, a, and outer ra- 
dius, b. Consider a position vector, x. In the original 
coordinate system (Fig.lA) it has components, x l , and 
in the transformed coordinate system (Fig. IB), x % . Of 




FIG. 1: The thick blue line shows the path of the same ray 
in (A) the original Cartesian space, and under two differ- 
ent interpretations of the electromagnetic equations, (B) the 
topological interpretation and (C) the materials interpreta- 
tion. The position vector x is shown in both the original and 
transformed spaces, and the length of the vector where the 
transformed components are interpreted as Cartesian compo- 
nents is shown in (C). 



course, its magnitude, r, is independent of coordinate 
system 
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where g^ji is the metric of the transformed space. In 
the materials interpretation, (Fig.lC), we consider the 
components, x % , to be the components of a Cartesian 
vector, and its magnitude is found using the appropriate 
flat space metric 
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Perhaps the simplest spherical cloak transformation 
maps points from a radius, r, to a radius, r', according 
to the following linear function 
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which we apply over the domain, < r < b, (or equiva- 
lently, a < r' < b). Outside this domain we assume an 
identity transformation. (All equations in the remainder 
of this article apply only to the transformation domain.) 
We must always limit the transformation to apply only 
over a finite region of space if we wish to implement it 
with materials of finite extent. Note that when r = 
then r' = a, so that the origin is mapped out to a fi- 
nite radius, opening up a whole in space. Note also that 
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when r = b then r' = b, so that space at the outer bound- 
ary of the transformation is undistortcd and there is no 
discontinuity with the space outside the transformation 
domain. 

Now since our transformation is radially symmetric, 
the unit vectors in materials interpretation and in the 
original space must be equal. 



where we have eliminated any dependence on the compo- 
nents of x in the original space, x l , or the magnitude, r. 
We can now drop the primes for aesthetic reasons, and we 
need not make the distinction between vectors and one- 
forms as we consider this to be a material specification in 
flat, Cartesian, three-space, where such distinctions are 
not necessary. Writing this expression in direct notation 
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Expressing the components of the position vector in the 
transformed space in terms of only the components in the 
original space, using (12), we obtain. 
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Now that we have this expression, we need not worry 
about the interpretations of transformed space, we can 
just proceed in standard fashion to compute the transfor- 
mation matrix. To take the derivative of this expression 
we note that 
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and obtain the transformation matrix 

ax l x k 5\6kj 



The components of this expression written out are 
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To calculate the determinant of this matrix we note that 
we can always rotate into a coordinate system such that 

(x i ) = (r,0,0) (18) 

then the determinant is, by inspection, given by 
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If we assume that our original medium is free space, then 
the permittivity and permeability will be equal to each 
other. As a short hand, we define a tensor, n l J , to 
represent both. 
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Though this definition is suggestive of refractive index, 
n % 3 would only represent the scalar refractive index 
if the permittivity and permeability were additionally 
isotropic, which is not the case here. 

Working out the algebra, we find that the material 
properties are then given by 
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where r <g> r is the outer product of the position vector 
with itself, also referred to as a dyad formed from the po- 
sition vector. We note, for later use, that the determinant 
can be easily calculated, as above, using an appropriately 
chosen rotation 

det^fAVY— V (23) 



B. Cylindrical Cloak 



To analyze a cylindrical cloak we will use two projec- 
tion operators. One which projects onto the cylinder's 
axis, (which we will call the third coordinate or 2-axis), 
and one that projects onto the plane normal to the cylin- 
der's axis. 
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We do not mean to imply that these are tensors. We de- 
fine these operators to perform these projections onto the 
third coordinate and the plane normal to the third co- 
ordinate in whatever basis (including mixed bases) they 
are applied to. Thus we will refer to their components 
with indices up or down, primed or un-primed, at will. 
We now use the transverse projection operator to define 
a transverse coordinate. 
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The coordinate transformation for the cylindrical case is 
the same as that of the spherical case in the two dimen- 
sions normal to the cylinder's axis. Along the cylinders 
axis the transformation is the identity. Thus we have for 
the transformation matrix. 
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or written in index form 



A 1 ' 



ap z p k 6\6 k j 



(27) 



4 



Again, we can easily calculate the determinant by rotat- 
ing into a coordinate system where 



(x*) = (?) = (p,o,o) 

then we find the determinant to be 
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The material properties in direct notation and dropping 
the primes are 
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Again we note the determinant for later use, which takes 
the rather simple form 
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III. HAMILTONIAN AND RAY EQUATIONS 

The Hamiltonian we will use for generating the ray 
paths is essentially the plane wave dispersion relation[17]. 
We derive it here, briefly, to show our choice of di- 
mensionality for the relevant variables. We begin with 
Maxwell's curl equations in SI units 
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We assume plane wave solutions with slowly varying co- 
efficients, appropriate for the geometric limit 

E = E e i(fe ° k - X_u ' t) H = J_H e i(fcok - x - wt) (33) 



Here ijo = y/Uo/^o is the impedance of free space, giv- 
ing E and H the same units, and k — w/c making k 
dimensionless. We use constitutive relations with dimen- 
sionless tensors e and /x. 
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Plugging (33) and (34) into the curl equations (32) we 
obtain 

kxE - /xH = k x H + eE = (35) 

Eliminating the magnetic field we find 

kx(/i- 1 (kxE o ))+eE o = (36) 

Defining the operator, K[13], 

Kik = ^ijkkj (37) 

the dispersion relation (36) can be expressed as a single 
operator on E , 



(K/x _1 K + e) E = 



(38) 



which must be singular for non-zero field solutions. The 
dispersion relation expresses that this operator must have 
zero determinate. 
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Now for material properties derived from transforming 
free space, e and ji are the same symmetric tensor, which 
we call n. In this case the dispersion relation has an 
alternate expression 
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This can be proved by brute force, evaluating the two ex- 
pression for a matrix, n, with arbitrary symmetric com- 
ponents, or perhaps some other clever way. The latter 
expression is clearly fourth order in k, but has only two 
unique solutions. Thus we discover that in media with 
s = fi the ordinary ray and extraordinary ray (found in 
anisotropic dielectrics) are degenerate. This can also be 
seen by noting that in free space the ordinary ray and 
extraordinary ray follow the same path. A coordinate 
transformation cannot separate these two paths, so they 
will follow the same path in the transformed coordinate 
space and thus also in the equivalent media. 

The Hamiltonian then easily factors into two terms 
that represent degenerate modes. Further it is easy to 
show ( by plugging (41) into (42) ) that the Hamiltonian 
may be multiplied by an arbitrary function of the spatial 
coordinates without changing the paths obtained from 
the equations of motion, (only the parameterization is 
changed), thus we can drop the factor, l/det(n), and 
our Hamiltonian is 



H = f (x) (knk - det (n)) 



(41) 



where / (x) is some arbitrary function of position. The 
equations of motion are[17] 
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where r parameterizes the paths. This pair of coupled, 
first order, ordinary differential equations can be inte- 
grated using a standard solver, such as Mathematica's 
NDSolve. 



IV. REFRACTION 



The equations of motion, (42), govern the path of the 
ray throughout the continuously varying inhomogenous 
media. At discontinuities, such as the outer boundary of 
a cloak, we must perform boundary matching to deter- 
mine the discontinuous change in direction of the ray, i.e. 
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refraction. Given ki on one side of the boundary we find 
k 2 on the other side as follows. The transverse compo- 
nent of the wave vector is conserved across the boundary. 



(ki 



o 



(43) 



where here n is the unit normal to the boundary. This 
vector equation represents just two equations. The third 
is obtained by requiring the wave vector to satisfy the 
plane wave dispersion relation of the mode represented 
by the Hamiltonian. 



H(k 2 ) = 



(44) 



These three equations determine the three unknowns of 
the vector components of k 2 . Since H is quadratic in 
k, there will be two solutions, one that carries energy 
into medium 2, the desired solution, and one that carries 
energy out. The path of the ray, dn/dr, determines the 
direction of energy flow, so the Hamiltonian can be used 
to determine which is the desired solution. The desired 
solution satisfies 

-•n>0 (45) 

if n is the normal pointing into medium 2. 

These equations apply equally well to refraction into 
or out of transformation media. Refracting out into free 
space is much easier since the Hamiltonian of free space 
is just, H = k • k — 1. 



V. CLOAK HAMILTONIANS 

We now show specific examples of ray tracing. Below 
we will choose a specific form for the Hamiltonian, plug 
in the material properties and display the derivatives of 
the Hamiltonian, for both the spherical and cylindrical 
cloak. 




FIG. 2: Rays traversing a spherical cloak. The transformation 
media that comprises the cloak lies between the two spheres. 




A. Spherical Cloak 



FIG. 3: Rays traversing a cylindrical cloak at an oblique an- 
gle. The transformation media that comprises the cloak lies 
in an annular region between the cylinders. 



For the spherical cloak (Fig. 2), the Hamiltonian which 
yields the simplest equations is 

lb 
2~ 
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Plugging in the material properties from (22) and (23) 
we obtain 
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Taking the derivatives, (which is straight forward partic- 
ularly in index form), yields 
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B. Cylindrical Cloak 

For the cylindrical cloak (Fig. 3), the Hamiltonian 
which yields the simplest equations is 
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Plugging in the material properties form (30) and (31) 
we obtain 



H = -kTk--^-^(p.k) 2 + - 



b(p-a) 



p(b-a) 
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(49) 

For taking the derivatives we note that the derivative of 
the transverse position vector with respect to the position 
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vector is the transverse projection operator. 
dx 



The derivatives are thus 
dH 2ap - a 2 , , . 
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these properties to perform ray tracing. Examples, of 
spherical and cylindrical cloaks are worked out in some 
detail. Some of the value in this effort is to provide in- 
dependent confirmation that the material properties cal- 
culated from the transformation do indeed cause electro- 
magnetic waves to behave in the desired and predicted 
manner. Eventually, this technique will become more ac- 
Zk cepted and independent confirmation will not be needed. 

One can see what the waves will do much more easily by 
applying the transformation to the rays or fields in the 
original space where the behavior is simpler if not trivial. 
However, one may still want to perform ray tracing on 
these media to see the effects of perturbations from the 
ideal material specification. 



(51a) 
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VI. CONCLUSION 



We have shown how to calculate the material proper- 
ties associated with a coordinate transformation and use 
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